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1.  Introduction 


The  fracture  and  failure  of  brittle  materials  has  remained  a  core  focus  in  Army  related  sciences. 
Of  particular  interest  is  the  destabilization  of  an  axisymmetric  solution  in  problems  that  are 
completely  axisymmetric.  This  destabilization  problem  largely  motivated  Grinfeld  and  Wright^  to 
study  a  model  for  dynamic  failure  that  depends  on  the  interplay  between  two  physical  effects: 
internal  elastic  energy,  and  the  energy  associated  with  breaking  chemical  bonds.  Because  of  these 
two  considerations,  the  model  is  referred  to  as  a  mechanochemical  model  for  fracture.  The 
motivation  of  the  mechanochemical  model  is  based  on  Stress  Driven  Rearrangements  Instabilities 
of  phase  interfaces  Grinfeld^  and  is  summarized  in  Grinfeld^  and  Kassner  et  al.."^  Motivated  by 
the  Stress  Driven  Rearrangements  Instabilities,  the  model  was  then  used  in  a  different  context  to 
form  the  mechanochemical  model,  Grinfeld  and  Wright,^  which  was  numerically  implemented  in 
MATLAB  to  handle  quasi-static  loading  cases.  Grinfeld  et  al.^  used  this  implementation  in 
MATLAB  to  established  the  appearance  of  radially  damaged  zones  and  the  destabilization  of  the 
axisymmetric  solution. 

A  clear  and  condensed  presentation  of  damage  theory  is  given  by  Kachanov,^  as  well  by 
Chaboche,’’  *  and  a  short  introduction  to  damage  as  an  internal  state  variable  can  also  be  found  in 
the  book  by  Holzapfel.^  In  the  types  of  damage  models  considered  by  Kachanov,  the  energy 
density  equation  is  modified  by  a  reduction  factor  or  a  damage  function  that  depends  on  an 
internal  state  variable,  the  damage.  As  one  might  expect,  the  damage  reduces  the  internal  energy 
thereby  reducing  the  stress  response  of  the  damaged  material.  In  the  mechanochemical  model 
discussed  here,  the  damage  also  contributes  to  the  internal  energy.  This  addition  to  the  energy 
potential  creates  an  interplay  between  mechanical  and  chemical  constituents  and  produces  the 
source  of  the  instability. 

The  original  intent  was  to  implement  the  mechanochemical  model  discussed  in  Grinfeld  and 
Wright^  into  the  government  code  SIERRA.  The  model  discussed  in  Grinfeld  and  Wright^  is 
based  on  a  modified  linear  elasticity.  In  this  report,  however,  we  focus  on  a  finite  deformation, 
nonlinear  formulation  of  the  problem.  This  change  was  required  by  constraints  introduced  within 
the  finite  element  solver  SIERRA.  In  SIERRA,  constitutive  model  designs  are  split  between  rate 
integrated  and  hyperelastic  formulations.  However,  the  features  of  the  mechanochemical  model 
make  the  decision  between  the  two  nontrivial.  The  mechanochemical  model  largely  depends  on 
the  internal  elastic  energy,  which  makes  it  a  good  candidate  for  the  hyperelastic  approach.  The 
desired  linear  elastic  response  on  the  other  hand,  would  be  simplest  to  implement  in  the  rate 
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integrated  formulation.  At  first,  a  rate  integrated  version  of  the  mechanochemical  model  was 
implemented  in  SIERRA.  This  initial  implementation  determined  the  linear  elastic  energy 
through  numerical  integration  but  was  prone  to  error.  As  a  compromise,  an  alternate  formulation 
of  the  model  was  developed  that  mimicked  the  linear  elastic  response  using  the  log-strain  tensor. 
This  method  benefited  from  an  accurate  calculation  of  elastic  energy  and  has  proved  to  be  more 
stable.  The  model  and  its  implementation  in  SIERRA  are  discussed  in  this  report.  Unlike  the 
previous  implementation  in  MATEAB,^  which  could  only  handle  quasi-static  loading,  this 
implementation  allows  for  the  dynamic  problem  of  failure  to  be  studied  on  a  massively  parallel 
computational  architecture. 

The  report  is  organized  as  follows.  We  discuss  the  general  physical  features  and  introduce  the 
system  of  equations  that  describe  the  damageable  material  in  section  2.  This  is  done  in  the 
framework  of  the  linear  elastic  theory  and  largely  recapitulates  the  work  of  Grinfeld  and  Wright.^*’ 
In  section  3  we  introduce  the  log-strain  tensor  and  derive  a  new  formulation  of  the 
mechanochemical  model  for  finite  deformations.  Section  4  presents  the  numerical  details  of  how 
the  damage  evolution  is  handled  and  discusses  some  of  its  subtleties.  In  section  5  we  verify  the 
implementation  of  our  model  in  SIERRA  using  some  single  element  tests.  A  few  example 
problems  that  illustrate  applications  of  the  mechanochemical  model  are  discussed  in  section  6, 
and  possible  future  alterations  to  the  model  are  discussed  in  section  7.  We  make  concluding 
remarks  in  section  8  and  derive  the  equations  of  damage  evolution  in  appendix  A.  Appendix  B 
includes  the  code  used  to  implement  the  model  in  SIERRA. 
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2.  Mechanochemical  Model:  Linear  Elastic  Framework 


In  this  section  we  largely  recapitulate  the  work  by  Grinfeld  and  Wright^*’  but  with  additional 
comments  to  motivate  the  following  section  where  we  introduce  a  finite  deformation 
mechanochemical  model.  We  begin  by  introducing  the  terms  considered  in  our  energy  density 
potential  from  which  we  derive  the  elastic  stress  and  the  equations  that  govern  the  chemical 
kinetics. 

We  take  the  total  energy  density  to  be  a  function  of  the  mechanical  deformation  (or  strain)  and  the 
damage.  Specifically,  we  assume  the  energy  can  be  decomposed  additively  into  a  mechanical  part 
and  a  chemical  part,  i.e., 

6  ^mechanical  “f  ^chemical  •  (1) 

Holding  the  damage  fixed,  the  derivative  of  the  energy  with  respect  to  the  strain  gives  the  stress. 


de 

d  strain 


=  stress . 


damage  constant 


Similarly,  holding  the  strain  constant  and  taking  the  derivative  of  energy  with  respect  to  the 
damage  gives  the  driving  force  for  damage  evolution. 


Strain  constant 


oc  —damage. 


Now  we  define  the  mechanical  and  chemical  terms  for  a  linear  elastic  material.  We  assume  the 
mechanical  energy  is  reduced  by  a  damage  function  0,  which  depends  on  the  damage  k  so  that  the 
mechanical  energy  is  given  by 

^mechanical  0(^)  ^elastic  •  (4) 

Here  Cmcchanicai  IS  the  mechanical  energy  of  a  damageable  material,  and  Ceiastic  is  the  elastic  energy 
associated  with  a  hypothetical  undamageable  material  undergoing  the  same  deformation.  We  also 
assume  the  damage  function  (p^n)  reduces  the  elastic  energy  with  increasing  damage.  Specifically, 
we  take  our  damage  function  0(/s;)  to  be 


n  ,  Cmin) 


l-(l-Cmin)^  ■■  K  <  K* 


:  K>  K* 


The  functional  form  of  equation  5  is  a  simple  linear  ramp  between  undamaged  and  a  maximally 
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damaged  material.  The  parameters  k*  is  the  maximal  damage  value.  If  while  updating  the  value  of 
K  the  damage  exceeds  k*,  it  is  simply  set  to  k*  (see  section  4).  The  parameter  Cmin  is  the  lower 
limit  to  the  reduction  factor  of  a  fully  damaged  material. 


Using  a  linear  elastic  theory  for  the  elastic  energy  term,  we  obtain  the  following  energy  density 
per  unit  reference  volume: 


^mechanical  ^ 


Here  Uij  is  the  linear  elastic  strain  tensor, 


1  /  dui  duj 
Uij  =  -  I  t:  h 


2  \  dxj  dxi  J 

where  Ui  is  the  displacement  vector  of  the  deformation. 


(6) 


(7) 


The  second  term  in  the  energy  density  is  the  chemical  energy, 
quadratic  and  of  the  form 


_  'C  /  o'.2 

^'chemical  2  '  ^  ' 


This  energy  is  assumed  to  be 


Here  is  a  chemical  constant  with  dimensions  of  stress  and  k,°  the  damage  associated  with 
minimum  energy. 


(8) 


Combining  these  two  terms  into  an  energy  density,  we  have  the  so-called  Kachanov-Lifshitz 
function. 

e{ui\j,  k)  =  n(t){K)  +  2  ~ 

Thus,  by  considering  changes  in  the  energy  density,  and  either  holding  the  strain  or  the  damage 
constant,  we  obtain: 


de 


dui 


ik 


l-2u 


UllSili:  T  Uik 


—  ^ik 


(10) 


and 


de 

Ok 


Uik  constant 


^^^^l^stic  “1“  "C  ^  ) 


dcj) 


2jj^  ■  1™^  U  U(^rn\n)U 


+  ^  (k  -  fi;°)  . 


(11) 

(12) 


The  stress  tensor  retains  its  usual  form,  i.e.,  it  is  still  defined  in  terms  of  the  derivative  of  the 
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energy  with  respect  to  the  strain  tensor;  however,  it  is  modified  by  the  damage  function 
Similarly,  a  generalized-stress  is  derived  from  the  chemical  energy.  It  has  units  of  stress  and 
drives  changes  in  the  damage.  This  generalized-stress  is  referred  to  in  the  literature  as  a  chemical 
potential.  The  following  two  equations  result  from  combining  the  stress  equation  10,  with 
momentum  balance,  and  the  chemical  potential  1 1  with  a  kinetic  reaction  law: 


P- 


d'^Uj  _ 

dt‘^  dxk  dxk 


2(I){k)p, 


^ll^ik  T  ^ik 


dt 


On' 


1  -2z/ 


V 


and 

+  ^  (k  -  K°) 


(13a) 

(13b) 


The  dimensions  of  K  are: 

[Time]  [Stress]  ’  ^ 

The  specific  choice  of  equations  13  contain  two  physically  questionable  features,  which  we 
discuss  in  later  sections.  Namely,  a  degrading  bulk  modulus  and  reversible  damage.  Reversible 
damage  used  in  this  report  is  not  consistent  with  the  work  by  Grinfeld  and  Wright^°  where  the 
derivative  of  the  damage  is  forced  to  remain  positive.  See  section  7.1  for  details  regarding  the 
minor  alteration  in  the  constitutive  model  that  disables  healing. 


3.  Mechanochemical  Model:  Finite  Deformation  Formulation 


As  discussed  in  the  introduciton,  the  original  goal  of  the  research  was  to  implement  the 
mechanochemical  model  for  brittle  fracture  described  in  the  previous  section  by  equations  5 
and  13  as  a  material  model  in  SIERRA.  However,  the  choice  to  implement  the  constitutive  model 
in  SIERRA  carried  with  it  a  number  of  constraints.  In  this  section  we  derive  a  constitutive  relation 
for  an  isotropic  hyperelastic  material  (hyperelastic  for  fixed  damage)  that  captures  the  same 
physical  features  described  in  the  previous  section.  This  does  not  represent  the  most  general 
nonlinear  model,  and  in  fact  is  a  special  case  that  was  chosen  to  mimic  the  linear  elastic 
mechanochemical  model  when  used  in  the  small  strain  limit. 

3.1  Finite  Deformations  and  the  Log-Strain  Tensor 

The  following  briefly  covers  some  modern  continuum  mechanics  and  introduces  the  log-strain 
tensor.  Two  good  references  on  modem  continuum  mechanics  are  the  text  books  by  Holzapfel,^ 
and  Truesdell  and  Noll.^^  The  Cauchy  stress  tensor  T  is  a  linear  transformation  of  a  direction 
vector  to  a  traction  force  in  the  current  configuration.  It  can  be  written  in  terms  of  the  spherical 
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part,  i.e.,  the  pressure  p,  and  its  deviatoric  part  T*: 

T  =  -pI  +  T*,  p=-\trT.  (15) 

o 

We  use  the  convention  that  normal  stress  components  are  positive  in  tension.  The  deviatoric  part 
of  a  tensor,  denoted  by  a  superscript  *,  is  defined  as  follows: 

A*  =  A-hv[A]I,  (16) 

where  I  is  the  identity  tensor.  The  deformation  gradient  F  is  a  linear  transformation  from  the 
reference  configuration  to  the  current  configuration,  which  can  be  expressed  in  terms  of  a 
properly  orthogonal  rotation  matrix  R  and  positive  definite  symmetric  matrices  U  or  V,  called 
the  right-  and  left-stretch  tensors,  respectively,  using  the  polar  decomposition 

F  =  RU  =  VR.  (17) 

Both  U  and  V  share  the  same  eigenvalues  A*  called  the  principal  stretches.  The  Jacobian  J  is  the 
determinant  of  the  deformation  gradient, 

J  =  detF  =  detV  ,  (18) 

and  is  equal  to  the  ratio  of  the  current  specific  volume  to  the  reference  specific  volume.  The  right 
and  left  Cauchy-Green  deformation  tensors,  respectively,  are  defined  as 

C  =  F^F  =  and 

B  =  FF'^  =  . 

Since  C  and  B  are  symmetric  and  positive  definite,  the  eigenvalues  of  C  and  B  are  Xf.  The 
principal  invariants  of  a  second-order  tensor  A  are  given  by 

I A  =  tr(A) ,  II A  =  tr(A2) ,  and  III  a  =  tr(A=^) .  (21) 

If  two  tensors  share  the  same  eigenvalues  one  can  equate  the  principal  invariants  of  one  to  the 
other.  Thus,  for  isotropic  elastic  materials  it  can  be  shown  that  the  energy  density  can  only  depend 
on  deformation  gradient  through  the  three  invariants  of  B  (or  equivalently  C).  In  this  manner  the 
energy  density  of  an  elastic  material  can  be  expressed  as  some  function, 

e  =  eB{lB,nB,inB).  (22) 


(19) 

(20) 
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There  are  many  choices  of  finite  deformation  strain  tensors  that  are  consistent  with  the  linear  elastic 
strain  tensor  in  the  infinitesimal  limit.  These  are  typically  expressed  in  terms  of  the  stretch  tensors, 
U  or  V,  or  the  Cauchy-Green  deformation  tensors,  C  or  B.  We  choose  to  use  the  log-strain  tensor, 
also  called  Hencky  strain.  This  particular  choice  of  strain  tensor  has  some  desirable  features  when 
compared  with  real  materials  at  moderate  strains and  is  often  considered  to  behave  similarly 
at  moderate  strains  to  an  infinitesimal  strain  measure.  The  log-strain  tensor  L  is  defined  as, 

L  =  \nV  =  ^\nB.  (23) 


Scheidler^^  considered  in  depth  the  log-strain  tensor  and  some  of  its  properties.  Notably  that 
tr  ly  =  In  J  and  hence  is  a  measure  of  the  volumetric  strain.  Furthermore,  the  deviatoric  part,  L* , 
of  L  is  independent  of  J  and  is  a  measure  of  shear  strain.  Since  II  =  trL  =  In  J  and  L*  is 
traceless,  the  energy  density  is  some  function  e'^  of  another  set  of  invariants: 


e  =  CLih,  IIl,  IIIl)  =  e'Jln  J,  • 


(24) 


Scheidler^^  using  the  log-strain  tensor,  also  showed  that  the  following  general  relationships  for  an 
isotropic  hyperelastic  material  can  be  derived  for  the  pressure  and  the  deviatoric  stresses: 


P  = 
JT*  = 


de{J,L*) 

de{J,  L*) 
dL*  , 


and 


(25a) 

(25b) 


The  log-strain  tensor  is  work  conjugate  to  the  Kirchoff  stress,  JT.  Work  conjugacy  follows  from 
a  general  kinematic  relationship,  from  which  it  is  shown  the  diagonal  components  of  the  time 
derivative  of  the  log-strain  tensor  are  equal  to  the  diagonal  components  of  the  rate  of  deformation 
tensor  in  the  principal  basis  for  the  left-stretch  tensor^^  and  the  assumption  of  an  isotropic  elastic 
material.  An  isotropic  elastic  implies  that  the  Cauchy  stress  is  also  diagonal  in  the  principal  basis 
for  the  left-stretch  tensor,  ensuring  that  the  inner  product  of  the  Cauchy  stress  and  the  rate  of 
deformation  tensor  is  equal  to  the  inner  product  of  the  Cauchy  stress  and  the  time  derivative  of  the 
log-strain  tensor. 


Only  the  second  and  third  principal  invariants  of  L*  are  nonzero,  so  that  the  deviatoric  stress  can 
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be  written: 


T 


[{L*fY. 


(26) 


* 


2  /  de 

1  {jih. 


L*  + 


3 

J 


The  deviatoric  stresses  depend  linearly  on  the  deviatoric  strains  through  the  second  invariant,  but 
their  dependence  on  the  third  invariant  is  quadratic  in  the  deviatoric  strains. Thus,  in  the  small 
shear-strain  limit  and  when  the  shear  modulus  is  independent  of  density,  one  can  conclude: 


T*  ^  ‘^L* ,  and 

fj 

p  ~  p{J) . 


(27a) 

(27b) 


Assuming  that  the  pressure  depends  only  on  the  Jacobian  J,  Scheidler^^  showed  that  the  strain 
energy  can  be  decoupled  additively  into  a  function  of  J  only  and  a  function  of  shear  only. 


e  6voi('7)  T  ('isoi^L  ) . 


(28) 


Taking  both  of  these  assumptions — specifically,  that  the  pressure  depends  only  on  the  Jacobian, 
and  the  deviatoric  stress  depends  linearly  on  the  deviatoric  strains — is  an  enabling  step  in  the 
formulation  of  the  finite  deformation  analog  to  the  linear  elastic  mechanochemical  model.  This 
formulation  is  similar  to  the  linear  elastic  theory  since  we  have  effectively  eliminated  the 
pressure-shear  coupling  that  would  otherwise  exist. 

3.2  Mechanochemical  Constitutive  Model 

We  now  define  the  decoupled  energy  density  for  the  mechanochemical  model.  As  before,  we  take 
the  energy  density  per  reference  volume  and  additively  decouple  it  into  a  mechanical  and 
chemical  part: 

6  ^mechanical  T  Cchemical  )  (29) 

and  express  the  mechanical  energy  in  terms  of  a  damage  function  0  and  the  elastic  energy: 

^mechanical  0(^)  ^'elastic  •  (30) 

The  elastic  energy  is  now  given  by  equation  28  and  the  total  energy  is 

e  0(^)  [Cvol('7)  T  ^isoi^L  )]  -f-  echem(^)  •  (31) 
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We  assume  the  equation  of  state  evoi(^)  is  given  by: 


evoi(</)  =  J  —  J  +  1)  ,  (32) 

where  JC  is  the  bulk  modulus;  other  equations  of  state  eould  be  used  in  the  future^  The  energy 
associated  with  the  deviatoric  deformation  is  taken  to  be 


eiUL*)  =  (33) 

so  that  the  total  energy  density  function  is  given  by 

e  =  (/.(k)  {/C(Jln  J-  J  +  1)  +  /itr  [(X*)^]}  +  ^  {k  -  .  (34) 

Equation  34  is  the  analog  of  the  Kachanov-Lifshitz  function  from  equation  9  consistent  with  finite 
deformations.  Using  equation  25  and  assuming  that  the  shear  strains  are  small  gives 


p  =  — 0(fi;)/Cln  J,  and 

Thus,  the  total  Cauchy  stress  tensor  is 


T  =  0(/^) 


/Cln  JI  + 


2/i 

T 


(35) 

(36) 


(37) 


The  dynamical  system  analogous  to  equations  13  for  the  finite  deformation  model  is 

d'^u 

p-^  =  V  ■  T ,  and  (38a) 

=  ^  (eeiastic(  J,  L*))  (38b) 

Note  that  Ceiastic  again  refers  to  the  elastic  energy  associated  with  a  hypothetical  undamageable 
material,  recall  equation  30.  Again  there  has  been  no  restriction  on  the  sign  of  the  derivative  of  k 
so  that  the  damage  is  reversible. 

^^Typically  the  (Jin  J  —  J  +  1)  term  in  the  expression  for  the  energy  that  relates  to  the  volumetric  contribution 
does  not  include  the  “+1”.  In  fact,  in  a  hyperelastic  model,  the  functional  form  of  the  energy  is  what  is  most  important 
within  an  arbitrary  constant  since  the  stresses  depend  on  derivatives  of  the  energy.  However,  since  we  now  would  like 
to  relate  this  energy  to  a  damage  growth  mechanism,  we  must  ensure  that  there  is  no  damage  growth  in  the  undeformed 
state  from  a  non-zero  energy  (i.e.,  setting  J  =  1  must  give  a  zero  energy). 
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4.  Implementing  the  Damage  Evolution  Equations 


In  this  section  we  discuss  how  the  coupled  dynamic  system  given  by  equation  38  is  handled  in  the 
code  and  how  it  is  different  from  the  case  where  the  left  hand  sides  of  equation  38  are  identically 
zero,  where  the  damage  develops  instantaneously.  Implementation  of  the  model  is  in  SIERRA 
SolidMechanics  (Presto/ Adagio  4.28;  Sandia  National  Laboratory)  a  Lagrangian  finite  element 
solver.  SIERRA  contains  both  an  explicit  solver  Presto  and  an  implicit,  nonlinear  preconditioned 
conjugant  gradient  solver.  Adagio.  The  implicit  code.  Adagio,  can  run  into  some  difficulties  since 
the  solution  can  become  nonunique  as  the  damage  increases.  In  other  words,  there  can  be  multiple 
deformations  that  correspond  to  the  same  stress,  and  the  solver  will  no  longer  converge.  Thus, 
while  the  code  works  with  Adagio,  its  use  is  not  recommended.  Since  SIERRA  is  a  solid 
mechanics  solver,  the  damage  model  must  be  introduced  without  disrupting  the  numerical 
stability.  Therefore,  there  are  some  additional  assumptions  discussed  here  that  are  needed  to  solve 
the  damage  evolution  alongside  the  mechanics. 

4.1  Dynamic  Evolution  of  Damage 

As  mentioned  in  section  2,  there  are  dimensions  of  time  in  the  kinetic  reaction  law  parameter  K. 
Thus,  there  are  two  time  scales  that  need  to  be  considered  when  solving  the  equations  of  motion. 
The  first  time  scale  is  related  to  wave  propagation  speeds  and  the  characteristic  dimensions  of  the 
body  and  the  second  is  how  quickly  damage  accumulates.  The  stability  of  the  explicit  code 
depends  critically  on  the  elastic  wave  speeds  and  the  minimum  element  size,  which  dictates  the 
time  step  used  during  the  calculations.  This  places  some  additional  requirements  on  the 
calculation  of  the  kinetic  reaction  law  since  it  too  needs  to  be  updated  at  each  time  step.  If  the 
kinetic  reaction  law  is  also  implemented  in  an  explicit  manner,  there  would  be  a  second  critical 
time  step.  Since  the  solutions  to  the  kinetic  reaction  law  are  exponentials  (see  appendix  A),  it 
could  become  an  incredibly  stiff  system  to  solve,  i.e.,  sensitive  to  error.  We  therefore  solve  the 
kinetic  reaction  law  for  an  incremental  change  in  time  so  that  we  can  directly  solve  for  the 
updated  damage  at  each  time  step. 

We  find  that  during  a  small  time  interval  /—)■/  +  At,  the  updated  damage  K{t  +  At)  can  be 
related  to  a  convolution  of  the  undamaged  elastic  energy  with  a  derivative  of  the  damage  function 
(see  appendix  A): 

nit  +  At)  =  K{t)e~^^^^  -  -  l)  -  K  |^(r)  .  (39) 
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We  assume  that  At  is  chosen  by  SIERRA  such  that  the  elastic  energy  is  varying  slowly.  This 
would  appear  to  be  a  safe  assumption  since  the  time  step  is  chosen  for  stability.  Thus,  we  pull  the 
elastic  energy  out  of  the  integral.  If  the  damage  function  is  given  by  equation  5,  we  can  also  pull 
the  derivative  of  the  damage  function  out  of  the  integral  and  evaluate  the  result  to  arrive  at  the 
following: 

k(«  +  At)  «  -  k”  -  l)  +  i  (l  -  ,  (40) 

s  ^ 

The  updated  damage  is  then: 


(A:  ( J  In  J  -  J  +  1)  +  /i  tr 

(41) 

where  superscripts  denote  the  time  step.  This  updated  damage  can  then  be  fed  into  the  stress 
calculation: 


T 


(/>(«)  /Cln  J/  + 


(42) 


The  minor  alteration  to  the  code  to  handle  the  case  where  the  damage  is  irreversible  is  discussed 
in  section  7.1. 


4.2  Instantaneous  Damage 

A  special  case  of  damage  evolution  is  when  the  left  hand  side  of  equation  38b  is  zero.  This  can  be 
achieved  by  taking  the  limit  as  K  goes  to  infinity  in  equation  40,  or  by  directly  solving 
equation  38b, 

K  =  K°  +  ^^—^^{}C{JlnJ-J  +  l)+^tr[{L*f])  .  (43) 

4.3  Additional  Details  of  the  Implementation 

SIERRA  provides  the  left-stretch  tensor  V  and  rotation  tensor  R  from  the  polar  decomposition  of 
the  deformation  gradient  tensor  F  at  the  so-called  next  time  step  n  +  1.  Erom  this  updated 
configuration,  one  must  determine  the  Cauchy  stress  in  the  unrotated  (reference  configuration). 

As  discussed  in  section  3,  the  model  depends  on  the  log-strain,  which  requires  taking  the  natural 
logarithm  of  a  tensor.  We  use  built-in  functions  in  SIERRA  to  perform  the  spectral  decomposition 
of  V,  from  which  the  log  can  be  taken.  The  overall  algorithm  is  outlined  as  follows: 


1.  Calculate  J  and  In  J  from  . 

2.  Determine  eigenvalues/vectors  for 
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3.  Calculate  the  log  of  the  eigenvalues  and  rotate  log  from  its  eigenbasis  back  to  current 
basis. 

4.  Calculate  from  In 

5.  Update  the  current  undamaged  elastic  energy, 

6.  Solve  for  the  new  damage  value  restrict  its  range  to  0  <  k  <  k*  using  equation  41 
or  43  as  selected  by  the  user. 

7.  Determine  the  value  of  the  damage  function 

8.  Determine  from  equation  42. 

9.  “Unrotate”  to  reference  configuration  by  calculating  R^TR. 

10.  Repeat  steps  1-9  for  all  elements  in  the  material. 


5.  Model  Verification 


This  section  compares  theoretical  predictions  of  the  model  against  simulation  results  from  single 
element  tests.  The  three  tests  considered  here  are  a  compression  test  where  the  volume  is  changed 
and  the  damage  evolves  instantaneously  to  the  equilibrium  state  so  that  equation  43  is  satisfied. 
The  second  test  is  a  shear  test  where  the  volume  does  not  change  and  again  the  damage  is  always 
instantly  updated.  The  last  test  is  a  compression  test  where  the  damage  evolution  is  time 
dependent,  i.e.,  according  to  equation  41.  In  all  tests,  we  chose  /C  =  1  Pa,  and  /i  =  0.75  Pa. 

5.1  Instantaneous  Damage  Evolution,  Confined-Compressive  Deformation 


We  compare  the  results  of  the  simulation  on  a  single  element  for  a  confined  compression  test.  The 
physical  components  of  the  deformation  gradient  F,  and  the  log-strain  tensor  In  V  imposed  are 


a 

0 

0  ■ 

In  a 

0 

0  ■ 

0 

1 

0 

,  [lnF]  = 

0 

0 

0 

_  0 

0 

1  _ 

0 

0 

0  _ 

(44) 
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The  Jacobian  J  =  deti^  =  a,  so  that  the  deviatoric  part  of  the  log-strain  tensor  is 


X*; 


I  In  a  0  0 

0  —  I  In  a  0 

0  0  —  llna 


(45) 


From  J  and  L*  we  can  solve  for  the  damage,  i.e.,  solving  equation  43  using  the  deformation 
given  in  equation  44.  The  solution  to  the  damage  is 


K  =  mm 


1,  ^  (  }C{alna  -  a  +  1)  +  ^ ]  (c^m  -  1) 


(46) 


This  results  in  the  magnitude  of  the  xx-component  of  the  Cauchy  Stress  to  be  given  by 


T  = 

XX 


/Clna  + 


4/i  In  a 
3q; 


mm 


f/C(alna-a  +  l)  +  ^^^^ 


ifimin  1) 


( ^rn.in. 


(47) 


The  results  of  the  theoretical  predictions  are  shown  in  figure  1  as  solid  lines  and  the  symbols  are 
values  taken  from  simulation  results.  Panels  a  and  c  plot  the  damage  k  as  a  function  of  the 
strain  X*^,  and  panels  b  and  d  plot  the  xx-component  of  stress  as  a  function  of  the  strain  X*^..  In 
these  plots,  =  0.02  and  k*  =  1.  The  top  two  plots  explore  the  result  of  holding  =  0.1  and 
varying  Cmin  from  0.2  to  1,  while  the  bottom  two  plots  explore  the  result  of  holding  c^m  =  0.3 
and  varying  ^  so  that  =  0.2  to  0.1.  In  all  cases  the  simulation  results  agree  completely  with 
the  theoretical  predictions. 


1)  +  1 


The  stress-strain  curves  in  panels  b  and  d  show  the  key  feature  required  for  the  physical 
instability.  Curves  that  exhibit  a  zero  slope  in  the  stress  followed  by  a  decrease  in  stress  at 
increased  strain  are  candidates  for  exploring  the  physical  instability.  We  leave  this  largely  without 
discussion  in  this  work  since  it  will  be  developed  further  in  subsequent  work. 
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Fig.  1  Single  element  compression  test.  Comparison  of  the  theoretical  prediction  (solid  lines)  against  the 
simulation  results  (symbols).  Panels  a  and  b  hold  fj,  =  0.1  and  vary  Cmin  from  0.2  to  1.  Panels  c  and  d 
hold  Cmin  =  0.3  and  vary  from  0.2  to  1.  Stress  made  positive  for  viewing  purposes  although  its  value  is 
negative  in  compression.. 
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5.2  Instantaneous  Damage  Evolution,  Isochoric  Shear  Deformation 


Here  we  eompare  the  results  of  the  simulation  on  a  single  element  for  an  isoehorie  shear 
deformation  in  the  Y-direetion  with  instantaneous  damage  evolution.  The  physieal  eomponents  of 
the  deformation  gradient  F  and  left-streteh  tensor  V  are 


■  1 

0 

0  ■ 

■  1 

0 

0 

0 

1 

a 

,  [^]  = 

0 

2+a^ 

a 

v^4+a^ 

_  0 

0 

1  _ 

_  0 

a 

2 

V  4+a2 

V 4+a^  -1 

(48) 


Using  MuPad  (Math Works),  one  ean  ealeulate  the  deviatorie  log-strain  tensor  L*,  the  damage  k, 
and  the  stress  T  from  the  previous  expressions  for  F  and  V.  Sinee  the  explieit  forms  of  these 
variables  are  quite  eomplieated,  they  are  not  reprodueed  here.  Instead,  figure  2  eompares  the 
simulation  results  (symbols)  against  the  theoretieal  predietion  (solid  lines).  This  figure  is  very 
similar  to  the  previous  one  with  the  exeeption  that  the  eomponent  of  the  stress  and  strain  that  is 
plotted  is  the  |/2;-eomponent.  The  figure  elearly  shows  agreement  between  theory  and  simulation 
and  that  the  physieal  instability  ean  manifest  itself  in  shear  as  well.^ 


^^The  level  of  shear  in  figure  2  likely  exceeds  the  range  of  validity  for  equation  27a. 
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Fig.  2  Single  element  shear  test.  Comparison  of  the  theoretical  prediction  (solid  lines)  against  the  simulation 
results  (symbols).  Panels  a  and  b  hold  =  0.1  and  vary  Cmin  from  0.2  to  1.  Panels  c  and  d  hold  Cmm  = 
0.3  and  vary  ^//r  from  0.2  to  1.. 
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5.3  Dynamic  Damage  Evolution,  Confined-Compressive  Deformation 


The  last  verification  step  presented  here  considers  the  time-dependent  damage  evolution.  We 
impose  a  very  rapid  compression,  identical  to  that  in  section  5.1  except  that  the  damage  is  now 
evolved  dynamically  according  to  equation  41.  Here  the  initial  value  of  the  damage  is  equal  to  the 
equilibrium  damage,  i.e.,  k(0)  =  kP  =  0.02,  and  the  kinetic  coefficient  K  =  0.0001  (1/Pa-s). 


For  the  following  analytical  calculation,  we  ignore  the  initial  transient  behavior.  We  assume  the 
body  is  initially  deformed  and  that  the  damage  starts  off  at  the  damage  associated  with  minimum 
energy  fi;(0)  =  Under  these  conditions  we  find  the  solution  to  the  damage  evolution  to  be 


K{t)  =  min 


1,  fi;(0)e 


-Kti 


K 


-Kti 


1) 


/r  {CfYiin  1) 

e 


(l  -  e  (  -/iloga  +  /C  (a  log 

(49) 


a 


with  a  corresponding  magnitude  of  the  xx  component  of  the  stress  evolution  given  by: 


q;  +  1) 


T  = 

XX 


K,  log  a  + 


4/r  log  a 
3a 


<PW))- 


(50) 


Figure  3  shows  good  agreement  between  the  simulation  and  the  theoretical  prediction.  The  only 
discrepancies  arise  from  the  initial  transient  behavior  that  we  ignored  (see  the  symbols  near  the 
vertical  axis  panels  b  and  d). 
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Fig.  3  Single  element  dynamic  compression  test.  Comparison  of  the  theoretical  prediction  (solid  lines) 
against  the  simulation  results  (symbols).  Panels  a  and  b  hold  =  0.1  and  vary  Cmin  from  0.2  to  1.  Pan¬ 
els  c  and  d  hold  Cmin  =  0.3  and  vary  ^  from  0.2  to  1.. 
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6.  Example  Problems 


The  previous  seetions  introdueed  the  equations  of  motion  and  the  verifieations  of  the  model’s 
implementation.  Here,  we  present  some  examples  of  how  the  model  eould  be  used. 

6.1  One-Dimensional  Dynamic  Confined- Compression 

In  this  seetion  we  eonsider  a  one-dimensional  (1-D)  problem  of  a  linear  alignment  of  material  that 
is  slowly  (eompared  to  the  elastie  wave  speeds)  eompressed  until  it  beeomes  unstable.  This 
problem  is  investigated  by  fixing  one  end  of  the  material  dX  Z  =  L  while  the  end  at  Z  =  0  is 
slowly  displaeed  aeeording  to  the  funetion: 

<5(t)  =  j  .  (51) 

1  ^max  .  t  ^ 

This  displaeement  funetion  has  a  diseontinuity  in  its  derivative  at  0  and  at  tc-  The  dimensions  and 
material  parameters  used  for  this  simulation  are  presented  in  table  1. 


Table  1  Summary  of  material  parameters  used  in  the  TD  stability  example.. 


L  (m)  p  (kg/m^) 

/C(Pa) 

/u(Pa)  aPa) 

K  (1/Pa-s)  Cmin 

tc  (S) 

^max  (m) 

1  le6 

10 

7.5  1 

0.1  0.3  0.02  1 

15000 

0.29 

Figure  4  shows  the  results  of  this  simulation.  Each  panel  of  the  figure  plots  the  damage 
distribution  versus  location  for  a  given  time.  The  difference  in  time  between  each  panel  is  not  the 
same.  At  first,  in  panels  a-e,  the  damage  evolution  is  similar  to  that  of  what  we  would  expect  in 
the  quasi-static  case,  i.e.,  it  is  proportional  to  the  strain  applied  to  the  system  and  shows  no 
localization  of  damage.  However,  in  panels/  and  g  the  damage  has  started  to  accumulate  at  the 
end  that  is  displaced.  This  localization  of  damage  culminates  by  panels  h  and  i  where  the 
saturation  of  damage  to  k*  in  panel  i  sends  off  a  relaxation  wave  in  the  damage.  This  wave 
propagates  toward  the  fixed  end  (Z  =  L)  where  a  reflection  of  the  wave  at  the  boundary  causes 
two  additional  localizations  of  damage  evolution,  panels  j  and  k.  At  this  point  the  damage 
evolution  slows  down  as  the  system  reaches  a  new  stable  equilibrium,  panels  l-o. 


19 


1 

0.8 

0.6 

0.4 

0.2 

0 


12350  s 


0.5  1 


1 

0.8 

0.6 

0.4 

0.2 

0 


13100  s 


0.5  1 


1 

0.8 

0.6 

0.4 

0.2 

0 


13230  s 


0.5  1 


Location  (m)  Location  (m) 


m)  14220  s 


1 . . i-i 

0.8 

0.6 _ L 

0.4 

0.2 

0 - - 

0  0.5  1 

Location  (m) 


1 

0.8 

0.6 

0.4 

0.2 

0 


17330  s 


0.5  1 

Location  (m) 


o)  24820  s 


1 . . 1—1'  ■ 

0.8 

0.6  _ 

0.4 

0.2 

0 - - 

0  0.5  1 

Location  (m) 


Fig.  4  A  One-Dimensional  alignment  is  slowly  compressed  from  its  Z  =  0  end.  Each  panel  plots  the 
damage  vs.  position  within  the  alignment,  where  the  plots  are  on  the  initial  configuration.  The  time  for  each 
panel  is  shown  in  the  upper  right  comer  and  the  time  difference  between  each  panel  is  not  constant.. 


One  subtle  feature  of  the  eonstitutive  model  is  that  the  ehemieal  energy  term  introduees  a 
ehemieal  potential  that  drives  the  damage  toward  This  means  that  the  material  ean  effeetively 
heal  and  that  the  damage  is  reversible.  This  feature  is  present  in  this  example  and  ean  be  seen  by 
eomparing  panels  e  and  o.  In  panel  e,  the  material  is  fairly  uniformly  damaged  with  a  value 
slightly  larger  than  0.6.  In  panel  o,  the  material  is  split  into  fully  damaged  material  k  =  1  and 
material  that  is  only  partially  damaged  k  ce  0.5.  In  some  eases  this  self-healing  feature  may  be 
undesirable,  and  a  modifieation  to  the  model  is  diseussed  in  seetion  7.1. 

6.2  Two-Dimensional  Radial  Damage  Bands 

In  this  seetion  we  eonsider  a  two-dimensional  (2-D)  problem  of  an  annular  dise  that  is  slowly 
(eompared  to  the  elastie  wave  speeds)  eompressed  until  it  beeomes  unstable  and  forms  radial 
damage  bands.  This  is  similar  to  the  problem  eonsidered  by  Grinfeld  et  al.,^  exeept  that  the 
meehanieal  system  is  not  evolved  quasi-statieally.  The  outer  radius  i?o  =  1  of  an  annular  dise  is 
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slowly  displaced  radially  toward  the  center  of  the  disc  according  to  the  piecewise  function  given 
in  equation  51.  The  inner  boundary  at  i?*  =  0.1  is  stress  free  and  the  disc  is  fixed  in  the  axial 
coordinate.  This  type  of  deformation  concentrates  stresses  near  the  inner  radius.  The  simulation 
used  the  material  parameters  given  in  table  2.  Because  of  the  choices  of  material  parameters  and 
the  imposed  deformation,  the  maximal  displacement  of  the  inner  radius  is  only  0.004  m  towards 
the  center. 


Table  2  Summary  of  material  parameters  used  in  the  2-D  radial  damage  band  instability  example.. 


Ri  (m) 

Ro  p  (kg/m^) 

/C  (kPa) 

M(kPa)  aPa) 

K  (1/Pa-s)  Cmin  H* 

(s)  ^max  (ni) 

0.1 

1  le6 

10 

7.5  1 

0.1  0.3  0.02  1 

3000  0.05 

Each  panel  in  figure  5  shows  the  damage  k  in  false-color  plotted  on  the  initial  configuration 
throughout  time.  The  time  elapsed  between  each  panel  is  175  s.  Initially,  the  damage  uniformly 
accumulates  throughout  the  material.  This  can  be  seen  by  noting  the  uniform  color  distribution  in 
panels  a-f.  The  solution  to  this  problem  without  damage  is  rotationally  symmetric  and 
concentrates  energy  at  the  small  opening.  This  gives  insight  into  why  the  damage  begins  to  localize 
near  the  center  and  where  the  material  instability  is  first  manifested.  The  accumulation  of  damage 
causes  a  local  weakening  near  the  opening  of  the  disc  and  the  germination  of  radial  damage  bands 
(panels  g  and  h)l.  Six  radial  damage  bands  clearly  form  (panel  i)  leaving  the  surrounding  material 
intact.  The  local  weakening  around  the  damage  bands  causes  a  strain  concentration  ahead  of  the 
existing  band.  This  subsequently  causes  the  instability  to  propagate  and  causes  the  band  to  grow 
predominantly  in  the  radial  direction  as  can  be  seen  in  panels  h-m.  The  length  of  the  band  grows  at 
a  finite  velocity,  where  one  can  estimate  the  speed  of  the  tip  of  the  vertical  (90°)  band  in  panels  h 
and  /,  0.001  m/s.  The  damage  is  arrested  as  the  stress  is  alleviated  and  the  system  finds  a  new  stable 
configuration.  No  additional  damage  accumulates  in  panels  m-o.  The  damage  bands  alleviate  hoop 
stresses  that  are  concentrated  at  the  center  opening  and  that  develop  from  the  imposed  deformation. 
In  future  developments  it  may  be  necessary  to  introduce  a  mechanism  for  hoop  damage  bands  to 
develop  from  radial  instabilities.  Two  possible  alterations  are  discussed  in  sections  7.1  and  7.2. 


Tn  these  damage  bands  the  material  is  still  intact.  In  an  abstract  way  these  bands  can  be  thought  of  as  a  proxy  for 
radial  cracks. 
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Fig.  5  Two-Dimensional  instability  and  the  development  of  radial  damage  bands.  A  2-D  annular  disc  is 
slowly  compressed  from  its  outside  edge  r  =  Rq.  Each  panel  plots  the  damage  in  false-color  as  it  varies 
with  position  within  the  disc.  Plots  are  on  the  initial  configuration.  The  elapsed  time  between  each  panel  is 
175  s.  At  first  the  damage  accumulates  symmetrically  before  a  material  instability  produces  the  localization 
of  damage  along  radial  lines.. 
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7.  Future  Considerations 


In  the  previous  sections  we  introduced  and  verified  the  implementation  of  a  mechanochemical 
model  for  failure.  In  this  section  we  consider  some  possible  alternatives  to  the  choices  of  the 
damage  function  0  and  the  energy. 

7.1  Removing  Healing  From  the  Mechanochemical  Model 

In  section  6.1  we  commented  that  the  quadratic  term  in  the  chemical  energy  can  result  in  a 
healing  process.  This  healing  process  is  driven  by  the  chemical  potential  and  may  not  be 
physically  realistic.  Analogous  to  Grinfeld  and  Wright/°  we  now  restrict  the  sign  of  k  so  that 
damage  is  irreversible.  Here  we  consider  some  ways  to  modify  the  current  model  and 
implementation  to  correct  for  this  feature. 

The  kinetic  evolution  of  the  damage  is  given  by  the  rate  equation  38b.  Repeated  here  for  clarity: 

=  ^  (eeiastic(7,  L*))  (52) 

For  our  specific  choice  of  0  (and  likely  in  more  general  cases),  the  derivative  of  the  damage 
function  with  respect  to  the  damage  is  negative.  Since  Ceiastic  >  0,  k  will  be  negative  if 

(53) 

J-  Cmin 

This  would  certainly  become  an  issue  in  a  test  where  the  loading  is  applied  very  rapidly  and  there 
are  elastic  waves  propagating  back  and  forth.  A  wave  propagating  past  would  initially  cause 
damage,  but  as  the  wave  continued  on,  the  damaged  material  might  be  in  a  lower  elastic  energy 
state  and  begin  to  self-heal  (recall  the  results  from  section  6.1). 

We  can  quite  trivially  remove  the  self-healing  feature.  The  algorithm  outlined  in  section  4 
calculates  the  mechanical  energy  associated  with  the  deformation.  This  energy  could  then  used  to 
evaluate  the  inequality  in  equation  53.  If  the  inequality  holds  true,  then  applying  the  traditional 
update  equation  41  would  result  in  a  smaller  damage  value  in  the  next  time  step.  However,  instead 
of  using  the  update  equation,  one  could  set  dn/dt  =  0,  and 

.  (54) 
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This  type  of  modification  might  be  sufficient  to  observe  radiai  and  hoop  damage  bands  a  dynamic 
test  where  heaiing  wouid  have  been  an  issue  before. 

7.2  Decoupling  Bulk  Damage  and  Deviatoric  Damage  Functions 

Here,  we  consider  a  simple  extension  to  the  current  model  where  the  only  modification  is  that 
there  are  two  damage  functions  that  can  affect  the  volumetric  and  isochoric  responses  separately: 

e  01  (^)Cvol('-^)  T  02(^)6iso(-^  )  T  6chem(^)  •  (55) 


The  stress  is  no  longer  given  by  equation  42,  but  instead 


2/r 


T  =  0i(fi;)/Cln  JI  +  . 

fj 


(56) 


The  damage  still  evolves  in  time;  however,  derivatives  of  0i  and  02  with  respect  to  k  will  enter 
into  the  kinetic  equation  separately: 


1  Ok  _  901  ^  ,  902 

Kdt~  dn  +  dK  ^  ^  ^  • 


(57) 


This  would  be  particularly  useful  if  the  material  was  more  sensitive  to  damage  in  shear  than  in  its 
bulk  response.  One  possible  choice  of  damage  functions  is 


01  (^)  1  (1  ^min)  r  and 

K 


02  (k)  =  max 


0,1  -  a  — 

K* 


(58) 

(59) 


A  dimensionless  parameter  a  >  0  sets  how  rapidly  the  damage  function  02  drops  to  zero.*  The 
maximum  operator  is  needed  here  since  the  original  notion  of  k*  is  lost  by  introducing  a  which 
would  otherwise  make  02  negative  when  an  >  k*.  In  this  example  Cmin  is  dropped  from  02  so 
that  the  material  completely  loses  its  deviatoric  response  and  behaves  more  like  a  soft,  nearly 
incompressible  material  when  it  is  maximally  damaged. 


An  alternate  approach  to  this  type  of  modification  would  be  to  introduce  a  second  damage 
parameter  K2,  with  its  own  damage  evolution  and  maximal  damage  In  this  case  there  would  be 
two  internal  state  variables  and  two  damage  evolution  equations  ki  to  solve.  This  also  would 
require  the  notion  that  damage  carries  a  direction  instead  of  simply  being  a  scalar. 

*  There  is  still  no  restriction  on  the  sign  of  the  time  derivative  of  the  damage.  Thus,  bulk  and  deviatoric  damage  can 
still  heal  and  that  a  merely  adjusts  the  rate  of  damage  accumulation  to  be  preferentially  bulk  or  preferentially  shear. 
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7.3  Damage  Growth  Asymmetry  in  Compression  and  Tension 

We  introduce  a  pathway  to  add  an  asymmetric  growth  of  the  damage  in  tension  and  compression. 
This  type  of  model  might  be  useful  to  ensure  that  a  material  does  not  lose  strength  in  compression 
and  to  promote  damage  growth  in  tension.  The  modified  energy  density  is 

e  =  0(k)  [atH{J  -  1)  +  acH{l  -  J)]  evoi(T)  +  0(k)  6chem(^  ) ,  (60) 


where  H  is  a  Heaviside  function: 


H{x)  = 


1  :  a;  >  0 
0  :  a;  <  0 


(61) 


The  construction  in  the  square  brackets  in  equation  60  acts  as  a  switch  between  the  two  energy 
scale  factors  at  and  ac.  H{J  —  1)  takes  on  the  value  1  when  the  material  is  in  tension,  and  if 
at  >  etc,  the  energy  is  increased  by  at  —  etc-  Assuming  this  form  of  the  energy  density  gives  the 
desired  behavior  in  the  kinetic  law 


la,H(J  -  1)  +  a,H(l  -  J)]  +  «  (k  -  k”) 


(62) 


The  deviatoric  stress  is 


JT*  = 


de 


dr 


=  2fiL*, 


and  the  pressure  can  be  obtained  from  applying  equation  25a  to  equation  60, 


(63) 


P  = 


de 


=  -^  {[atH{J  -  1)  +  acH{l  -  J)]  /C  (Jin  J  -  J  +  1))  . 


(64) 


Evaluating  the  derivatives  and  simplifying  gives: 


p  =  [atH{J  -  1)  +  a^il  -  J)]  /C  In  J , 


(65) 


Now  the  slope  of  the  pressure  has  a  discontinuity  as  we  cross  J  =  1  from  above  or  below: 
dv  1 

^  =  [atH{J  -  1)  +  acH{l  -  J)]  /C-  +  [at  5{J  -  1)  -  etc  5(1  -  J)]  /C  In  J ,  (66) 

CJ  fj  fJ 


where 


lim  ^=a,}C 
j^i-  dJ 


(67) 
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and 


(68) 


lim  ^  =  at  1C. 
j^i+  oJ 

This  is  only  a  minor  concern  since  the  derivative  of  the  stress  response  already  has  a  diseontinuity 
from  the  damage  funetion  (see  figures  1  and  2). 


8.  Conclusions 


In  eonelusion  we  have  developed  a  meehanoehemieal  model  for  dynamie  failure  using  a  finite 
deformation  formulation  analogous  to  the  model  developed  by  Grinfeld  and  Wright.^  We  have 
eondueted  single  element  tests  to  verify  the  model  is  working.  We  have  also  outlined  some 
example  problems  that  the  model  was  implemented  to  study.  We  also  diseussed  some  of  the 
limitations  and  possible  future  alterations  that  eould  be  made  to  the  model. 
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In  the  main  text,  we  derived  a  kinetic  law  relating  the  evolution  of  the  damage  to  the  elastic 
energy  of  the  material, 

%  =  (^eelastic  +  e  («  -  •  (A-1) 

For  the  special  case  discussed  in  the  main  text,  equation  A-1  readily  integrates.  This  appendix 
considers  a  more  general  approach  that  could  be  used  for  other  damage  functions  0. 


For  notational  simplicity  we  write, 

and  rewrite  equation  A- 1  as 


dK 

dt 


-K  {g{n(t))  eeiastic(f)  +  ^  («  -  «°))  • 


(A-2) 


(A-3) 


Recall  that  the  elastic  energy  is  the  hypothetical  elastic  energy  of  an  undamageable  material,  i.e., 
the  energy  that  only  depends  on  the  displacements  and  the  elastic  moduli.  In  general,  g  might  not 
be  a  simple  function  of  the  damage  so  a  method  of  Laplace  transforms  is  useful  when  attempting 
to  solve  this  differential  equation.  In  the  following,  C  [f{t)]  =  F{s)  denotes  the  Laplace 
transform  of  f{t).  Similarly,  C~^  [-^(-s)]  =  f{t)  denotes  the  inverse  Laplace  transform  of  F{s). 


Applying  the  Laplace  transform  to  both  sides  of  equation  A-3  gives 


sC  [K{t)]  -  k(0)  =  -KC  [g{K{t))  eelastic(f)]  -  K^C  [K{t)]  + 


(A-4) 


which  is  an  algebraic  equation  that  can  be  solved  for  C 


/tfO)  1  / 


s  +  Ki  Ki\  s  s  +  Ki 
Applying  the  inverse  Laplace  transform  gives: 

«(0) 


KC  [g{n{t))  eelastic(f)] 
s  +  Ki 


K{t)  =  C  ^ 


K{t)  =  K{d)e  —  tF  (e 


c-^ 


K 


c 


-1 


K 


s  +  Ki 


C 


-1 


KC  [g{K{t))  eelastic(f)] 
s  +  Ki 


l)  -  KC  ^  [C  [g{K{t))  eeiastic(f)]  c  [e  . 


(A-5) 


(A-6) 

(A-7) 


The  last  term  is  the  inverse  Laplace  transform  of  the  product  of  two  Laplace  transforms,  which  is 
the  convolution  of  the  two  original  functions: 


K(t)  =  K(0)e  -  K°  (e  -  l)  -  iF  [  eeiastic(r)e  . 


(A-8) 
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Since,  in  general,  dcjy/dn  will  depend  on  k,  this  is  an  integral  equation  for  nit).  It  can  be  used  to 
update  the  damage  parameter  for  any  damage  function  as  a  convolution  of  84)/ dn  and  Ceiastic- 
Considering  a  small  time  step  At,  however,  one  can  derive: 

K{t  +  At)  =  -l)  -  K  eeiastic('r)e"^^^‘+'^*“^^dr  .  (A-9) 

If  At  is  chosen  for  stability  in  the  explicit  codes,  the  change  in  elastic  energy  will  likely  be  small 
over  the  time  interval  t  — )■  t  +  At.  Thus,  we  take  the  elastic  energy  outside  of  the  convolution 
integral.*  We  also  assume  that  the  damage  is  not  changing  rapidly  during  this  time  interval  so  that 
84) /dn  is  slowly  varying  and  can  also  be  taken  outside  the  integral.  In  the  special  case  of  4  given 
by  equation  5  this  holds  trivially,  since  84)/ 8k  is  constant.  Thus,  we  can  write  the  incremental 
change  in  the  damage  as: 


K{t  +  At)  ^  -  K-  -  1)  -  i  (1  -  Celastic  (A-10) 

Thus,  the  updated  damage  K{t  +  At)  depends  on  the  current  damage  K{t)  the  equilibrium  damage 
K°  and  the  elastic  energy. 

Using  a  superscript  n  denotes  that  the  variable  is  evaluated  at  the  current  time  step,  and  n  +  1  the 
future  time  step  gives 

-  K°  -  1)  -  ^  (1  -  •  (A-1 1) 


When  0  is  given  by  equation  5, 


84)  _  1  Cfnin 

8k  K* 


(A-12) 


and 

-  K°  -  l)  +  ^  (l  -  •  (A-13) 

4  K* 

In  the  text,  the  elastic  energy  term  is  used  in  two  different  ways.  First,  that  of  a  linear  elastic 
material  and  later,  a  hyper-elastic  material. 


*In  future  models  where  d(j)/dK  is  not  constant,  additional  state  variables  may  be  needed  in  the  constitutive  model, 
which  keeps  track  of  past  elastic  energies  and  past  damage  values.  These  might  be  used  in  some  Simpson  quadrature 
rule  to  simplify  the  convolution  so  that  it  could  be  evaluated  numerically. 
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^include  <models  /  B  r i tt  1  e E 1  a s  t i  c  .  h> 

^include  <cmath> 

# include  <kinematics/lame_general_methods  .h> 
^include  <kinematic s  /  LameUtil  .  h> 

^include  <kinematics/ Kinematics  .h> 
using  namespace  std; 
namespace  lame  { 


Material  *  BrittleElastic : : createMaterial(  const  MatProps  &  props  )  { 
return  new  BrittleElastic (  props  ) ; 

) 

//  There  are  three  state  variables  (although  only  one  changes  right  now) 

//  The  first  two  are  basic,  BULK.MODULUS  and  SHEAR_MODULUS 
//  The  third  which  evolves  in  time  is  the  damage  parameter  KAPPA 
//  To  make  this  model  work,  however,  you  need  to  specify 

//  XI_CHEMICAL  damage  chemical  constant  ,  K_KINETIC ,  KAPPA_ZERO  the  steady  state  damage 
//  Also  will  need  C_MIN  and  KAPPA_STAR  which  appear  in  the  damage  function  phi  . 

//  These  are  all  constant  throughout  the  run  and  are  properties. 

//  This  is  the  constructor  for  the  BrittleElastic  model.  The  properties 
//  it  reads  in  are  the  BULK_MODULUS  and  SHEAR_MODULUS .  These 
//  are  stored  in  the  properties  array  . 

// 

BrittleElastic :: BrittleElastic( const  MatProps  &  props)  : 

Material (props )  { 
setFlag(USE_LEFT_STRETCH) ; 
setHyper ( ) ; 

// 

//  Material  Property  Definitions 
//  What  is  read  in  from  the  input  dec 

//  YOUNGS_MODULUS,  POISSONS_RATIO ,  XI_CHEMICAL,  K_KINETIC ,  KAPPA_ZERO,  C_MIN,  KAPPA_STAR,  -e- 
INSTANT  —  8  total  properties 


mat_name  =  "BRITTLEELASTIC"  ; 
num_material_properties  =  8; 
initializeProperties () ; 

setMaterialProperty (0  ,  " YOUNGS JVIODULUS"  , props)  ; 
setMaterialProperty ( 1  ,  "POISSONS_RATIO"  , props)  ; 

//  The  default  behavior  is  linear  elastic  without  damage 
setMaterialPropertyDef ault ( 2 ,0.0) ; 
setMaterialProperty (2  ,  " XI_CHEMICAL "  , props)  ; 
setMaterialPropertyDefault (3  ,0.0) ; 
setMaterialProperty (3 , "K_KINETIC" , props) ; 
setMaterialPropertyDefault (4  ,0.0) ; 
setMaterialProperty (4  ,  "KAPPA_ZERO"  , props )  ; 
setMaterialPropertyDefault (5  ,0.0) ; 
setMaterialProperty (5 , "C_MIN" , props) ; 
setMaterialPropertyDefault (6 ,1.0) ; 
setMaterialProperty (6  ,  "KAPPA_STAR"  , props)  ; 
setMaterialPropertyDefault (7 ,0) ; 
setMaterialProperty (7  ,  "INSTANTANEOUS"  , props)  ; 
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LAME_FORTRAN(elastic_property_check) (  num_material_properties , 
&properties [0]  ); 


// 

//  State  Variable  Definitions 

//  We  are  using  5  state  variables.  The  bulk  and  shear  are  two.  then  damage  kappa,  then  two 
energies 

// 

num_state_vars  =  5; 

set_state_variable_alias  ( "YOUNGS_MODULUS"  ,  0); 
set_state_variable_alias  ( "POISSONS_RATIO"  ,  1); 
set_state_variable_alias  ( "KAPPA_INTERNAL"  ,  2)  ; 
set_state_variable_alias  ( "MECHANICAL_ENERGY"  ,  3)  ; 
set_state_variable_alias  ( "CHEMICAL_ENERGY"  ,  4); 


int  BrittleElastic :: initialize (  matParams  *  p  )  { 
const  double  ym  =  properties [0] ; 
const  double  nu  =  properties  [  1  ] ; 

originalYoungs  =  properties [ 0 ] ; 
originalNu  =  properties  11]; 


xi  =  properties  1 2 ] ; 
bigk  =  properties  1 3  ] ; 
kzero  =  properties [4] ; 
cmin  =  properties [ 5 ] ; 
kstar  =  properties [6] ; 
instant  =  properties  1 7 ] ; 


const  double  kappa  =  properties  1 4  ] ;  //initially  we  are  neutral  /  steady  state  damaged 

double  stateOld  =  p— >state_old ; 

double  ^  stateNew  =  p— >state_new ; 

for  (int  i(0);  i  <  p— >nelements ;  ++i)  { 

// 

//  Set  the  default  state  to  the  block  elastic  constants 

// 


StateOld [0] 
stateOld [ I ] 
state01d[2] 
stateOld [ 3 ] 
stateOld [4] 
stateOld  += 
StateNew  += 


=  stateNewtO] 

=  stateNewtl] 

=  stateNewt2] 

=  stateNewl3] 

=  stateNew[4] 

5; 

5;//5  internal 


ym ; 
nu ; 

kappa ; 

0.0; 

0.0; 

state  variables 


return  0; 

} 


// 

//  The  getStress  method  returns  the  stress  given  the  strain 
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//  and  the  time  step  . 

// 

int  BrittleElastic : : getStress (  matParams  *  p  )  { 
const  double  dt  =  p— >dt ; 

//  const  double  *  strainRate  =  p— >strain_rate  ; 
double  *  rotation  =  p— >rotation; 
const  double  *  leftStretch  =  p— >left_stretch ; 
const  double  *  stressOld  =  p— >  stress_old; 
const  double  *  stateOld  =  p— >  state_old; 
double  ^  stressNew  =  p— >  stress_new; 
double  *  stateNew  =  p— >  state_new; 

const  double  onehalf  =  1. 0/2.0; 
const  double  onethird  =  1. 0/3.0; 

double  devlnV  [6]; 
double  evals  [ 3 ] ; 
double  evecs  [ 9 ] ; 
double  LQT  19]; 
double  resultmat  [9]; 
double  tracelnV; 

for  (int  i(0);  i  <  p— >nelements ;  ++i)  { 

const  double  ym  =  originalYoungs ;  // stateOld  [0]  <  0  ?  0  ;  stateOld  [0];  //  ym  must  be  ^ 

positive 

const  double  pr  =  originalNu;  // stateOld  1 1  ]  <  —0.9999999  ?  —0.9999999  : 
const  double  bulk  =  ym/3.0/(  1 .0  — 2.0*pr)  ; 
const  double  shear  =  ym/ 2 . 0/ ( 1 . 0  +  pr )  ; 

const  double  kappa  =  state01d[2]  >  kstar  ?  kstar  : 

state01d[2]  <  0  ?  0  :  state01d[2];  //keep  kappa  less  than  kstar  and  greater  than  0 

const  double  J  =  lef  tStretch  [xx  ]*(  leftStretch  [yy]  *  left  Stretch  [  z  z]  — leftStretch  [  yz  ]* 
leftStretch  tyz ] ) 

— leftStretch  txy]  h=  ( leftStretch  txy]  *  leftStretch  [  zz]  — leftStretch  [  zx]  *  leftStretch  [yz  ]  ) 

+  le  ft  St  retch  1  zx]  =i=  ( leftStretch  txy  ]  *  lef  t  St  retch  [yz]  — leftStretch  [yy]  =!=  leftStretch  [  zx  ]  )  ; 
const  double  InJ  =  FastLog(J); 
const  double  sheartwoover J  =  shear*2.0/J; 
const  double  InJdivJ  =  InJ/J; 


devlnV[xx]= 
devlnV[yy]= 
devlnV [ zz ] = 
devlnV[xy]= 
devlnV[yz]= 
devlnV [ zx] = 


leftStretch [yy " 
leftStretch [ zz " 
leftStretch [xy " 
leftStretchlyz " 
leftStretch [ zx " 


LAME_FORTRAN ( lame_eigen) ( 1 ,1 , devlnV, evals , evecs ) 
//Then  use  these  values  to  calculate  logV  in  its 


;//calculate  eigen  values  and  vectors 
eigenbasis  . 


//then  convert  basis  back  to  where  we  are  using:  QXQ^  =  V 
evals [0]=FastLog(evals [0] ) ; 
evals [ 1 ] =FastLog(evals [ 1 ] )  ; 
evals [2] =FastLog(evals [ 2 ] ) ; 


36 


161 

162 

163 

164 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

180 

181 

182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 

201 

202 

203 

204 

205 

206 

207 

208 

209 

210 

211 

212 

213 


//make  Lambda*Q^T,  name  is  LQT  . 

//  this  is  taking  the  rows  of  Q^T  and 
LQT[fxx]  =  evals  [0]  evecs  [ fxx ] ;  // 
LQT[fxy]  =  evals [0] * evecs [ fyx ] ;  // fyx 
LQT[fxz]  =  evals  [0]  evecs  [ f zx ] ; 


mult  by 
because 


lambda_i 
it  is  Q^T, 


for  each 
not  Q 


row 


LQT[fyy] 

LQT[fyz] 

LQT[fyx] 


evals [ 1 ] *  evecs [ fyy ] ; 
evals [ 1 ] *  evecs [ f zy ] ; 
evals [ 1 ] *  evecs [ fxy ] ; 


LQT[fzz] 

LQT[fzx] 

LQT[fzy] 


evals [2]^ evecs [ f zz ] ; 
evals [2] *  evecs [ fxz ] ; 
evals [2]^ evecs [ fyz ] ; 


//now  Q  matrix  is  just  the  evecs  matrix,  so  multiply 
//  InV  =  evecs  *  LQT 

//this  lame  call  to  fortran  does  not  work,  so  I  replaced  it  with  a  direct  calculation  instead 
//  LAME_FORTRAN  ( lame_tensor_product3636){l  ,1  ,  evecs  ,LQT,  resultmat);/  / multiply  evecs  *LQT 
//result  mat  should  be  symmetric  since  it  is  now  In  V 


resultmat [ fxx] 
resultmat [ fxy] 
resultmat [ fxz ] 


evecs [ fxx] *  LQT [ f xx]  +  evecs [ fxy] *  LQT [ fyx]  + evecs [ f xz ] *  LQT [ f zx ] ; 
evecs [ fxx] *  LQT [ f xy ]  +  evecs [ fxy] *  LQT [ fyy ]  + evecs [ f xz ] *  LQT [ f zy ] ; 
evecs [ fxx] *  LQT [fxz ]  +  evecs [ fxy] *  LQT [ fyz ]  +  evecs [ f xz ] *  LQT [ f zz  ] ; 


resultmat [ fyx] 
resultmat [fyy] 
resultmat [ fyz ] 


evecs [ fyx] *  LQT [ f xx]  +  evecs [ fyy] *  LQT [ fyx ]  + evecs [ f yz ] *  LQT [ f zx ] ; 
evecs [ fyx] *  LQT [ f xy ]  +  evecs [ fyy] *  LQT [ fyy ]  + evecs [ f yz ] *  LQT [ f zy  ] ; 
evecs [ fyx] *  LQT [fxz ]  +  evecs [ fyy] *  LQT [ fyz ]  +  evecs [ f yz ] *  LQT [ f zz ] ; 


resultmat [ f zx] 
resultmat [ f zy ] 
resultmat [ f zz ] 


evecs [ f zx] *  LQT [ f xx]  +  evecs [ f zy ] *  LQT [ fyx ]  + evecs [ f zz ] *  LQT [ f zx ] ; 
evecs  [  f  zx]  *  LQT  [  f  xy  ]  +  evecs  [  f  zy]  LQT  [  fyy ]  + evecs  [  f  zz  ]  *  LQT  [  f  zy  ]  ; 
evecs [ f zx] *  LQT [fxz ]  +  evecs [ f zy ] *  LQT [ fyz ]  +  evecs [ f zz ] *  LQT [ f zz ] ; 


t race In V  =  resultmat [ fxx ]+ resultmat [ fyy ]+ resultmat [ f zz ] ; 


devlnV  [ XX ]  =  resultmat  [  f xx]— one t hi rd=)« trace InV ; 
devlnV  [yy]  =  resultmat  [  f yy]— one t hi rd=)« trace InV ; 
devlnV [  zz ]  =  resultmat  [  f  z z]—onet hi rd=)« trace InV ; 
devlnV [xy] = resultmat [ fxy ] ; 
devlnV [yz ]=resultmat [fyz ] ; 
devlnV [ zx] = resultmat [ f zx ] ; 


//  update  the  mechanical  energy 

stateNew[3]=bulk*(  J*lnJ—J+l  )+shear*(de  vlnV[xx]*  devlnV  [xx]  + devlnV  [yy]*devlnV[yy]  +  devlnV[zz]*-<-^ 
devlnV [ zz ] + 

2.0*  devlnV [xy] *  devlnV [xy]  +  2 .0*  devlnV [yz ] *  devlnV [yz ]  +  2 .0*  devlnV [ zx] *  devlnV [ zx ] ) ; 

//  update  the  new  damage 

const  double  expbxt  =  (instant  ==  1)  ?  0  :  exp(— bigk*xi*dt )  ; 

double  kappaNew  =  kappa  *  expbxt— k  zero  *  (expbxt  —  1.0)  —  (1 .0  — expbxt )  *  ( cmin—  1.0)  /kstar  *stateNew  [  3  ]  /-e^ 
xi ;  //  alphaNew  has  a  mu  in  it 

kappaNew  =  kappaNew  >  kstar  ?  kstar  : 

kappaNew  <  0.0  ?  0.0  :  kappaNew; 

const  double  phi  =  1.0  —  (1.0  —  cmin) *kappaNew/kstar ; 
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//  finally  calculate 
stressNew[0]  = 
stressNew[l]  = 
stressNew[2]  = 
stressNew[3]  = 
stressNew[4]  = 
stressNew[5]  = 


the  cauchy  stress 

phi5)=bulk*ln  J  +  phi^sheartwoover  j5i<devlnV[xx  ] ; 
phi5t=bulk*lnJ  +  phi^sheartwoover  J*devlnV[yy  ] ; 
phi*bulk*ln J  +  phi shear twooverJ*devlnV [  zz  ] ; 

phi  shear  t  woover  J*devlnV[  xy 
phi^f'sheartwoover  J*devlnV[  y  z 
phi=f«sheartwoover  J*devlnV[  zx 


/I  XX 

//yy 
//  zz 
;  //xy 

;  //  yz 

;  //  zx 


//rotate  (or  unrotate)  back  to  reference  RFtR 
//in  Sierra  T  is  called  sigma 

//use  resultmat  to  temporarily  store  the  info 
unRotate ( stressNew , rotation ) ; 


stateNew[0]  =  ym^phi ;// original  youngs  modulus  is  reduced  by  phi 
stateNew[i]  =  pr ; // poissons  ratio  is  not  affected  by  damage 
stateNew[2]  =  kappaNew; 
stateNew[3]  =  stateNew [3 ] *phi ; 

stateNew[4]  =  onehalf  *  (kappaNew  —  kzero)  >f'(kappaNew  —  kzero)*xi; 
rotation  +=9; 

stressOld  +=  6; 
stressNew  +=  6; 
leftStretch  +=6; 
stateOld  +=  5; 

StateNew  +=5; 


return  0; 


int  BrittleElastic : : getInitialElasticModuli (  matParams  *  p  )  { 
const  double  *  stateOld  =  p— >  state_old; 
for  (int  i(0);  i  <  p— >nelements ;  ++i)  { 

p->pcVars_old[PC_YOUNGS_MODULUS  ]  [  i]  =  stateOldlO];  //  9 .0  s  t  ate  0  Id  1 0]  =i=  s  t  ate  0  Id  [  1  ]  /  ( 3 . 0  * 
stateOldlOJ+stateOld  [1])  ; 

p— >pcVars_old[PC_POISSONS_RATIO]  [  i]  =  stateOldll];  // ( 3 . 0  =i=  s  t  ate  0  Id  [0]  —  2.0*  s  t  at  eO  Id  1 1  ] ) -^^ 
/{ 2.0 *(3.0*  StateOld  [0]  +  StateOld  [  1])  )  ; 

} 

return  0; 

) 

int  BrittleElastic :: pcElasticModuli (  matParams  *  p  ) 

{ 

const  double  *  stateOld  =  p— >  state_old; 
double  *  StateNew  =  p— >  state_new; 
double  *  oldE  =  p->pcVars_old [PC_YOUNGS_MODULUS ] ; 
double  *  newE  =  p->pcVars_new [PC_YOUNGS_MODaLUS ] ; 
double  *  oldNu  =  p->pcVars_old[PC_POISSONS_RATIO ] ; 
double  *  newNu  =  p->pcVars_new [PC_POISSONS_RATIO ] ; 

for  (int  i(0);  i  <  p— >nelements ;  ++i) 

{ 

oldE[i]  =  stateOld[0]; 
newE[i]  =  stateNew[0]; 
oldNu[i]  =  stateOldll]; 
newNu[i]  =  stateNew[l]; 


38 


267 


stateOld  +=  5; 

268  stateNew  +=5; 

269  } 

270  return  0; 

271  } 

272  void  BrittleElastic  ::  unRotate ( double  *  stressNew,  double  *  rotation) 

273  { 

274  double  resultmat  [6]; 

275  //  R'^SR 

276  resultmat  [xx]  =  rotation  [  f xx  ]  ^  ( rotation  [  f xx]  *  stressNew  [xx  ]  +  rotation  [  f yx]  *  stressNew  [  xy-e^ 

]  +  rotat  ion  [  fzx]  *  stressNew  [  zx  ] )  +  rotation  [  fyx  ]  ( rotation  [  fxx]  stressNew  [xy]  + 

rotation [ fyx] *  stressNew [yy]  + 

277  rotation  [  fzx]  *  stressNew  [  yz  ] )  +  rotation  [  fzx  ]*(  rotation  [  fyx]  *  stressNew  [yz  ]  +  rotation[fxx]>f'-<-^ 

stressNew [ zx ]  +  rotation [ fzx] * stressNew [ zz ] ) ; 

278  resultmat  [yy  ]  =  rotation  [  f xy  ]  *  ( rotation  [  f xy  ]  *  stressNew  [xx  ]  +  rotation  [  f yy  ]  *  stressNew  [  xy-e^ 

]  +  rotat  ion  [  fzy]  *  stressNew  [  zx  ] )  +  rotation  [  fyy  ]*(  rotation  [  fxy]  stressNew  [xy]  + 
rotation[fyy]>f'StressNew[yy]  + 

279  rotation  [  fzy]  *  stressNew  [yz  ] )  +  rotation  [  fzy  ]*(  rotation  [  fyy]  stressNew  [yz  ]  +  rotation[fxy]*-<-^ 

stressNew [ zx ]  +  rotation [ fzy] * stressNew [ zz ] ) ; 

280  resultmat  [  zz  ]  =  rotation  [  f xz  ]  *  ( rotation  [  f xz  ]  *  stressNew  [xx  ]  +  rotation  [  f y z  ]  *  stressNew  [  xy-<-^ 

]  +  rotat  ion  [  fzz  ]*  stressNew  [  zx  ] )  +  rotation  [  fyz  ]*(  rotation  [  fxz  ]=)«  stressNew  [xy]  + 
rotation  [  fyz  ]  stressNew  [yy]  + 

281  rotation  [  f  z  z  ]  *  stressNew  [  y  z  ] )  +  rotation  [  f  zz  ]  *  ( rotation  [fyz]^  stressNew  [yz  ]  +  rotation  [  f  xz  ]  *  •<-^ 

stressNew[ zx]  +  rotation [ fzz ]* stressNew [ zz ] ) ; 

282 

283  resultmat  [xy  ]  =  rotation  [  f xx  ]  *  ( rotation  [  f xy ]  *  stressNew  [xx  ]  +  rotation  [  f yy  ]  *  stressNew  [  xy-<-^ 

]  +  rotat  ion  [  fzy]  *  stressNew  [  zx  ] )  +  rotation  [  fyx  ]*(  rotation  [  fxy]  stressNew  [xy]  + 
rotation  [  fyy]  stressNew  [yy]  + 

284  rotation  [  f  zy  ]  *  stressNew  [  y  z  ] )  +  rotation  [  f  zx  ]  *  ( rotation  [  fyy  ]  *  stressNew  [yz  ]  +  rotation  [  fxy  ]  * 

stressNew[ zx]  +  rotation [ fzy] * stressNew [ zz ] ) ; 

285  resultmat  [yz  ]  =  rotation  [  fxy  ]  ^  ( rotation  [  f xz  ]  *  stressNew  [xx  ]  +  rotation  [  f y z  ]  *  stressNew  [  xy-<-^ 

]  +  rotat  ion  [  fzz  ]*  stressNew  [  zx  ] )  +  rotation  [  fyy  ]  ( rotation  [  fxz  ]=)«  stressNew  [xy]  + 

rotation [ fyz ]* stressNew [yy]  + 

286  rotation  [  f  z  z  ]  *  stressNew  [  y  z  ] )  +  rotation  [  f  zy  ]  *  ( rotation  [  fyz  ]  *  stressNew  [yz  ]  +  rotation  [  fxz  ] 

stressNew[ zx]  +  rotation [ fzz ]* stressNew [ zz ] ) ; 

287  resultmat  [  zx]  =  rotation  [  f xz  ]  >)«  ( rotation  [  f xx]  *  stressNew  [xx  ]  +  rotation  [  fyx]  *  stressNew  [  xy-e^ 

]  +  rotat  ion  [  fzx]  *  stressNew  [  zx  ] )  +  rotation  [  fyz  ]  sf' ( rotation  [  fxx]  stressNew  [xy]  + 
rotation [ fyx] *  stressNew [yy]  + 

288  rotation  [  f  zx]  *  stressNew  [  y  z  ] )  +  rotation  [  f  zz  ]  *  ( rotation  [  fyx]  *  stressNew  [yz  ]  +  rotation  [  f  xx]  * 

stressNew[ zx]  +  rotation [ fzx] * stressNew [ zz ] ) ; 

289 

290  for(unsigned  int  j=0;  j<6;++j) 

291  stressNew  [  j  ]  =  resultmat  [  j  ] ; 

292  } 

293  }  //  lame 
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